library(ggplot2)
require(data.table)
library("dplyr")
library("ggpubr") 
require(tidyverse) #for easy data manipulation and visualization
require(caret) #for easy machine learning workflow
library(ncdf4)
library(raster)
library("forcats")
library(RColorBrewer)

#change to full local path to the replication data
wd <- "C:\\Users\\YFan\\OneDrive - NORCE\\NFood_replication\\"

case = c("RCP45","RCP45_85LU")
casename = c("RCP45",expression('RCP45'[(SSP5LU)]))
region = c("NAM", "SEA", "EUA", "SAM", "CAM", "SAF", "Global")

#load intermediate data prepared by "Data.Fig.S5-S6.R"
load(paste0(wd,"/data/RData/Fig.S5-6.LUC-Analysis-WeightedQuantiles.Rdata"))


#=========Supplementary Fig.S5 Crop yield change (a), regional climate shift (b: air temperature, c: precipitation) 
#=========and change in soil fertility (d) for each crop from SSP5 to SSP2 (SSP2-SSP5) under ER
df.ssp2 <- df3.all[Region=="Global" & Year!=2100 & Scenario=="RCP45"]
df.ssp5 <- df3.all[Region=="Global" & Year!=2100 & Scenario=="RCP45_85LU"]
df <- df.ssp2[Variable=="YIELD"]

df$yield.dif=(df.ssp2[Variable=="YIELD"]$Value-df.ssp5[Variable=="YIELD"]$Value)/df.ssp5[Variable=="YIELD"]$Value*100
df$tsa.dif=df.ssp2[Variable=="TSA"]$Value-df.ssp5[Variable=="TSA"]$Value
df$ppt.dif=(df.ssp2[Variable=="PRECIP"]$Value-df.ssp5[Variable=="PRECIP"]$Value)/df.ssp5[Variable=="PRECIP"]$Value*100
df$sminn.dif=(df.ssp2[Variable=="SMINN"]$Value-df.ssp5[Variable=="SMINN"]$Value)/df.ssp5[Variable=="SMINN"]$Value*100
#df$sand.dif=(df.ssp5[Variable=="SAND"]$Value-df.ssp2[Variable=="SAND"]$Value)/df.ssp2[Variable=="SAND"]$Value*100
#df$clay.dif=(df.ssp5[Variable=="CLAY"]$Value-df.ssp2[Variable=="CLAY"]$Value)/df.ssp2[Variable=="CLAY"]$Value*100
df$fertn.dif=(df.ssp2[Variable=="FERTN"]$Value-df.ssp5[Variable=="FERTN"]$Value)/df.ssp5[Variable=="FERTN"]$Value*100
df$irrig.dif=df.ssp2[Variable=="IRRIG"]$Value-df.ssp5[Variable=="IRRIG"]$Value


jpeg(paste0(wd,"/figures/supplementary/Fig.S5a.jpeg"), units = "in", width = 6, height = 4, res = 600)

p <- ggplot(df, aes(Year, yield.dif, group=Crop)) + geom_line(aes(color = Crop)) +
  scale_x_continuous(breaks =c(2006, 2050, 2099)) +
  scale_y_continuous(expression(Delta~Yield*" (%)")) #(MtC year"^-1*")"), 
# Change the color palette
#p <- set_palette(p, palette = get_palette("Set5", 6))
##remove grid and background or remove all
p + theme_bw() +  #theme_classic() +  #theme_minimal()
  geom_hline(yintercept = 0) + # add horizontal reference line
  ggtitle("(a)") + theme(plot.title = element_text(hjust = 0.5)) +  
  theme(title = element_text(size = 14)) +
  theme(legend.title = element_blank(), legend.text = element_text(size=14)) + #, legend.position = "bottom") +
  theme(axis.title=element_text(size=14), axis.text = element_text(size=12), strip.text.x = element_text(size = 12)) +
  theme(axis.title.x =element_blank(), # axis.text.x = element_text(angle=45,hjust=1), 
        panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  theme(axis.title.y = element_text(color = "black"),
        axis.title.y.right = element_text(color = "blue"),
        axis.text.y.right = element_text(color = "blue"))
dev.off()


jpeg(paste0(wd,"/figures/supplementary/Fig.S5b.jpeg"), units = "in", width =6, height =4, res = 600)
p <- df %>% 
  ggplot(aes(Year, tsa.dif, group=Crop)) + geom_line(aes(color = Crop)) +
  #geom_line(aes(Year, ppt.dif/20, group=Crop, color=Crop), linetype=2) + #color=fct_reorder(Crop, effect, fun=median, .desc = TRUE))) +
  scale_x_continuous(breaks =c(2006, 2050, 2099)) +
  scale_y_continuous(expression(Delta~TSA*" (°C)")) #limits = c(-1,1), #breaks =c(-0.9, 0, 0.9), #all sec axis use scale from primary
#sec.axis = sec_axis(~ .*20, name = expression(Delta~PRECIP*" (%)"))) 
##remove grid and background or remove all
p + theme_bw() +  #theme_classic() +  #theme_minimal()
  geom_hline(yintercept = 0) + # add horizontal reference line
  ggtitle("(b)") + theme(plot.title = element_text(hjust = 0.5)) +  
  theme(title = element_text(size = 14)) +
  theme(legend.title = element_blank(), legend.text = element_text(size=14)) + #, legend.position = "bottom") +
  theme(axis.title=element_text(size=14), axis.text = element_text(size=12), strip.text.x = element_text(size = 12)) +
  theme(axis.title.x =element_blank(), # axis.text.x = element_text(angle=45,hjust=1), 
        panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  theme(axis.title.y = element_text(color = "black"),
        axis.title.y.right = element_text(color = "blue"),
        axis.text.y.right = element_text(color = "blue"))
dev.off()


library("forecast")
df = df[, ppt.dif.ma:=ma(ppt.dif,10), by = .(Crop)]
jpeg(paste0(wd,"/figures/supplementary/Fig.S5c.jpeg"), units = "in", width =6, height =4, res = 600)
p <- df %>% 
  ggplot(aes(Year, ppt.dif.ma, group=Crop)) + geom_line(aes(color = Crop)) +
  #geom_line(aes(Year, irrig.dif, group=Crop, color=Crop), linetype=2) + #color=fct_reorder(Crop, effect, fun=median, .desc = TRUE))) +
  scale_x_continuous(breaks =c(2006, 2050, 2099)) +
  scale_y_continuous(expression(Delta~PRECIP*" (%)")) #limits = c(-1,1), #breaks =c(-0.9, 0, 0.9), #all sec axis use scale from primary
#p <- set_palette(p, palette = get_palette("Set5", 6)) # Change the color palette
##remove grid and background or remove all
p + theme_bw() +  #theme_classic() +  #theme_minimal()
  geom_hline(yintercept = 0) + # add horizontal reference line
  ggtitle("(c)") + theme(plot.title = element_text(hjust = 0.5)) +  
  theme(title = element_text(size = 14)) +
  theme(legend.title = element_blank(), legend.text = element_text(size=14)) + #, legend.position = "bottom") +
  theme(axis.title=element_text(size=14), axis.text = element_text(size=12), strip.text.x = element_text(size = 12)) +
  theme(axis.title.x =element_blank(), # axis.text.x = element_text(angle=45,hjust=1), 
        panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  theme(axis.title.y = element_text(color = "black"),
        axis.title.y.right = element_text(color = "blue"),
        axis.text.y.right = element_text(color = "blue"))
dev.off()


jpeg(paste0(wd,"/figures/supplementary/Fig.S5d.jpeg"), units = "in", width =6.7, height =4, res = 600)

df = df[, fertn.dif.ma:=ma(fertn.dif,5), by = .(Crop)] 
p <- df %>% 
  ggplot(aes(Year, fertn.dif.ma, group=Crop)) + geom_line(aes(color = Crop)) +
  #ggplot(aes(Year, sminn.dif, group=Crop)) + geom_line(aes(color = Crop)) +
  geom_line(aes(Year, sminn.dif, group=Crop, color=Crop), linetype=2) + #color=fct_reorder(Crop, effect, fun=median, .desc = TRUE))) +
  scale_x_continuous(breaks =c(2006, 2050, 2099)) +
  scale_y_continuous(expression(Delta~FERTN*" (%)"), #limits = c(-1,1), #breaks =c(-0.9, 0, 0.9), #all sec axis use scale from primary
                     sec.axis = sec_axis(~ .*1, name = expression(Delta~SMINN*" (%)"))) #, breaks =c(0, 10))) 
##remove grid and background or remove all
p + theme_bw() +  #theme_classic() +  #theme_minimal()
  geom_hline(yintercept = 0) + # add horizontal reference line
  ggtitle("(d)") + theme(plot.title = element_text(hjust = 0.5)) +  
  theme(title = element_text(size = 14)) +
  theme(legend.title = element_blank(), legend.text = element_text(size=14)) + #, legend.position = "bottom") +
  theme(axis.title=element_text(size=14), axis.text = element_text(size=12), strip.text.x = element_text(size = 12)) +
  theme(axis.title.x =element_blank(), # axis.text.x = element_text(angle=45,hjust=1), 
        panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  theme(axis.title.y = element_text(color = "black"),
        axis.title.y.right = element_text(color = "blue"),
        axis.text.y.right = element_text(color = "blue"))
dev.off()


#=====global FERTN trend same as Fig.1e Nitrogen fertilization
# df.fertn <- df3.all[Variable=="FERTN" & Region=="Global",]
# df.fertn[,fertn:=Value]
# df.fertn[,Variable:=NULL]
# df.fertn$area <- df3.all[Variable=="AREA" & Region=="Global",]$Value
# df.fertn.gl <- df.fertn[,.(fertn.gl=weighted.mean(fertn, area)), by=.(Scenario,Year,Region)]
# p <- df.fertn.gl %>%
#   ggplot(aes(Year, fertn.gl, group=Scenario)) + geom_line(aes(color = Scenario)) +
#   scale_x_continuous(breaks =c(2025, 2050, 2075)) +
#   guides(color = guide_legend( keywidth = 1, keyheight=0.8, label.position = "right", label.hjust = 0, ncol=1))
# p + theme_bw() + theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank())+
#   ggtitle(expression('Land use')) + theme(plot.title = element_text(hjust = 0.5,size=13)) +
#   theme(legend.text = element_text(size=11), legend.title = element_blank(), legend.background = element_blank(),
#         legend.position = c(0.76, 0.77)) +
#   theme(axis.title=element_text(size=11), axis.text = element_text(size=11)) +
#   theme(panel.spacing = unit(0, "lines")) +  #reduce facet spacing
#   theme(strip.text = element_text(size = 11))




#=========Supplement Figure S6 for land use change analysis of ER scenario===========


###########function for new boxplot style
plot.wtQ <- function(df, index){
  print(paste("df is for dataframe of variable", var," ;", "index is the plot number", index))
  p0 <- df %>% 
    #mutate(Crop = factor(Crop, levels=rev(levs))) %>%
    mutate(Scenario = factor(Scenario, levels=case, labels=c("SSP2", "SSP5"))) %>%
    ggplot() + coord_flip()+ #horizontal bars
    geom_rect(aes(xmin=Xmin, xmax = Xmax, ymin = Q25, ymax = Q75, fill = Scenario), colour = NA)+
    geom_rect(aes(xmin=Xmin+0.3, xmax = Xmax-0.3, ymin = Q5, ymax = Q95, fill = Scenario), colour = NA)+
    geom_rect(aes(xmin=Xmin+0.4, xmax = Xmax-0.4, ymin = Q1, ymax = Q99, fill = Scenario), colour = NA)+
    geom_rect(aes(xmin=Xmin, xmax = Xmax, ymin = Q50, ymax = Q50, fill = Scenario), colour = "white", size=1)+
    facet_wrap(~Crop, nrow=3, labeller=label_parsed, scales='free_x', strip.position = "top") + #edit facet position and name
    #scale_fill_brewer(palette="Paired") + 
    scale_fill_manual(name="", values= alpha(c("#E69F00", "#56B4E9"),0.5))+ #light red and blue
    scale_x_discrete(labels=c("SSP2", "SSP5")) +
    scale_y_continuous(unit) 
  p0 <- p0 + theme_bw()+  #theme_classic() +  #theme_minimal()
    ggtitle(" ") + theme(title = element_text(size = 14)) + #ggtitle(paste(index, region[n]))
    guides(fill = guide_legend(label.position = "left", ncol=1,label.hjust = 0, keyheight=0.5, keywidth=2)) +
    theme(legend.title=element_blank(),legend.text = element_text(size=12), legend.position = "left") +
    theme(axis.title=element_text(size=14), axis.text = element_text(size=12), strip.text.x = element_text(size = 14)) +
    theme(#axis.title.x =element_blank(), # axis.text.x = element_text(angle=45,hjust=1), 
          panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
    theme(plot.margin = margin(t = 0, r = 0.5, b = 0, l = 0, unit = "cm"))
  return(p0)
}

#=======Fig.S6.weighted percentiles of regional yeild and climate
#select one crop per region with all variables as facets at y axis (modified facet style)
region
#crop.select <- c("Soy", "Rice",  "Wheat", "Sugarcane", "Maize", "Cotton") #pick one crop for "NAM" "SEA" "EUA" "SAM" "CAM" "SAF"
crop.select <- c("Maize","Rice","Wheat", "Soy", "Sugarcane", "Cotton")
# var.s <- c("YIELD",  "TSA", "PRECIP", "SMINN")
# varunit = c("'Yield (tC ha'^-1*')'", "TSA*' (°C)'", "'PRECIP (mm yr'^-1*')'", "'SMINN (gN m'^-2*')'")
var.s <- c("YIELD", "SMINN", "TSA", "RH")
varunit = c("'Yield (tonnes ha'^-1*')'", "'SMINN (gN m'^-2*')'", "TSA*' (°C)'", "'RH (%)'")

for (n in 1:6) { #one crop per region
  crop.s <- crop.select[n]
  df0 <-  df.all[Crop%in%crop.s & Region==region[n] & Variable%in%var.s]

  df0[Variable=="YIELD", WtQ1:=WtQ1*10^6/0.45*0.85] #MtC/ha to tonnes/ha
  #df0[Variable=="TSA", WtQ1:=sprintf("%.0f", WtQ1)] #no decimal for TSA (does not work)
  df1 <-  df0[Probs=="5%"]
  df1$Q1 <- df0[Probs=="1%"]$WtQ1
  df1$Q5 <- df0[Probs=="5%"]$WtQ1 
  df1$Q25 <-df0[Probs=="25%"]$WtQ1
  df1$Q50 <-df0[Probs=="50%"]$WtQ1
  df1$Q75 <-df0[Probs=="75%"]$WtQ1
  df1$Q95 <-df0[Probs=="95%"]$WtQ1
  df1$Q99 <-df0[Probs=="99%"]$WtQ1
  df1[Scenario=="RCP45_85LU", Xmin:=1]
  df1[Scenario=="RCP45_85LU", Xmax:=2]
  df1[Scenario=="RCP45"]$Xmin <- 2+1/2
  df1[Scenario=="RCP45"]$Xmax <- 3+1/2
  
  #==calc 2075-2099 mean yield of rainfed and irrigated crops in tonnes/ha 
  yield.irrig<-df3.all[Crop%in%crop.s & Variable=="YIELD.IRRIG" & Region==region[n] & Year>=2075, 
                       .(Variable="YIELD", MEAN.IRRIG=sum(Value)/25*10^6/0.45*0.85), by=.(Crop, Scenario)]
  yield.rain <-df3.all[Crop%in%crop.s & Variable=="YIELD.RAIN" & Region==region[n] & Year>=2075, 
                       .(Variable="YIELD", MEAN.RAIN=sum(Value)/25*10^6/0.45*0.85), by=.(Crop, Scenario)]
  df2 <- df1[yield.irrig, on=.(Crop, Scenario)]
  df2 <- df2[yield.rain, on=.(Crop, Scenario)]
  df2 <- df2[Variable!="YIELD", MEAN.IRRIG:=NA]
  df2 <- df2[Variable!="YIELD", MEAN.RAIN:=NA]
  
  case <- c("RCP45", "RCP45_85LU")
  #plot.wtQ(df1[Variable==var], "(a)") #MtC/ha to tC/ha
  p2 <- df2 %>% 
    mutate(Variable = factor(Variable, levels=var.s, labels=varunit)) %>%
    mutate(Scenario = factor(Scenario, levels=case, labels=c("SSP2", "SSP5"))) %>%
    ggplot() + coord_flip()+ #horizontal bars
    #geom_rect(aes(xmin=Xmin+0.4, xmax = Xmax-0.4, ymin = Q1, ymax = Q99, fill = Scenario), colour = NA)+
    geom_rect(aes(xmin=Xmin+0.4, xmax = Xmax-0.4, ymin = Q5, ymax = Q95, fill = Scenario), colour = NA)+
    geom_rect(aes(xmin=Xmin+0.1, xmax = Xmax-0.1, ymin = Q25, ymax = Q75, fill = Scenario), colour = NA)+
    geom_rect(aes(xmin=Xmin, xmax = Xmax, ymin = Q50, ymax = Q50, fill = Scenario), colour = "white", size=0.8)+
    geom_point(aes(x=(Xmin+Xmax)/2, y=MEAN.RAIN, color = Scenario), size=2, shape=16)+ 
    geom_point(aes(x=(Xmin+Xmax)/2, y=MEAN.IRRIG, color = Scenario), size=1., shape=2, stroke=1)+ #triangle for irrigated
    facet_wrap(~Variable, nrow=2, labeller= label_parsed, scales='free_x', strip.position = "bottom") + 
    #scale_fill_brewer(palette="Paired") +
    scale_fill_manual(name="", values=alpha(c("#E69F00", "#56B4E9"),0.8)) + #light/transparent red and blue
    scale_color_manual(name="", values=c("#D55E00", "#0072B2"), guide=FALSE) + #dark red and blue
    scale_x_discrete(labels=c("SSP2", "SSP5")) + scale_y_continuous("")
    #scale_y_continuous("", labels = scales::number_format(accuracy = 0.1)) #control decimal place fpr all variables
  
  p2 + theme_bw()+  #theme_classic() +  #theme_minimal()
    ggtitle(paste(region[n],":",crop.s)) + theme(title = element_text(size = 11), plot.title = element_text(hjust =0.5)) +
    guides(fill = guide_legend(label.position = "left", ncol=1,label.hjust = 0, keyheight=0.5, keywidth=2)) +
    theme(legend.title=element_blank(),legend.text = element_text(size=9), legend.position = "left") +
    theme(axis.title=element_text(size=10), axis.text = element_text(size=9), strip.text.x = element_text(size = 10)) +
    theme(axis.title.y =element_blank(), # axis.text.x = element_text(angle=45,hjust=1), 
          panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
    theme(plot.margin = margin(t = 0, r = 0.5, b = 0, l = 0, unit = "cm")) +
    theme(panel.spacing.x = unit(1, "lines")) +
    theme(strip.background = element_blank(), strip.placement = "outside")
  #together with strip.position = "left" or "bottom" to get traditional x/y axis title/units for facets
  
  ggsave(paste0(wd,"/figures/supplementary/Fig.S6.",region[n],"_",crop.s,".jpeg"), units = "in", width =4, height =3, dpi = 600)
  
}
 
 









